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This is a tutorial paper which highlights the 
influence of the computer not only on the modus 
operandi of circuit design but also on network 
theory itself. It reviews the topological proper- 
ties of linear graphs and describes a matrix- 
topological formulation of the network problem. 

In addition to the classical mesh, node, and cutset 
methods, a mixed method of analysis is described 
which is applicable to d-c, a-c, and transient 
problems. 


Numerical methods of solving linear and non- 
linear d-c network problems are discussed and 
a new approach to a~c analysis, using the mixed 
method and a numerical solution of the matrix 
eigenvalue problem, is described. The extension 
of this method to the transient analysis of linear 
networks is also explained. Finally, the problem 
of instability in the numerical integration of 
differential equations is discussed and several 
means of solving the problem are outlined, 


The Impact of Computers 


As pointed out in a review article by Kuo, 7 
and more recently by Christiansen,“ the computer 
is exerting a profound influence on the modus 
operandi of circuit design. With the aid of existing 
computer programs, the circuit designer may 
carry out ‘numerical experiments" which would be 
difficult or even impossible to duplicate with actual 
hardware on the bench. Often, he may be able to 
do such experiments faster and cheaper on the 
computer. He may analyze circuits for which cer- 
tain components have not yet been manufactured — 
or even invented— and he may vary parameters, 
either singly or in combinations, seeking for an 
optimum design. Moreover, he may do all this 
without the risk of burning up expensive or hard- 
to~get components even if he inadvertently over- 
drives the circuit with a megavolt input signal, 


Thus, the digital computer provides circuit 
designers with an analytical tool of great power 
and versatility. Generally speaking, the computer 
can be trusted to produce reliable solutions to an 
arbitrary network problem, except when the prob- 
lem involves pathological numerical difficulties 
or when the method of analysis has been poorly 


chosen, It is the user's responsibility, however, 
to pose a problem that accurately represents the 
circuit which the designer has in mind, since the 
computer does only what it is explicitly instructed 
to do and has no way of second-guessing the de- 
signer's real intent. It is obvious therefore that 
high-fidelity network models for diodes and tran- 
sistors will be needed, We shall not attempt to 
discuss the modelling of devices in this paper, 
but their importance cannot be ignored, 


In addition to its obvious impact on circuit 
design, the computer is also having a more subtle 
influence on network theory itself. For example, 
because matrix notation is so convenient in pro- 
gramming a digital computer, many existing net- 
work analysis programs are based on a matrix 
formulation of the network problem, This, together 
with the inherent elegance of the matrix approach, 
is helping to promote its acceptance among network 
theorists. 


In another direction, computers are helping 
to hasten the recognition of the inherently broad 
scope of network theory by drawing attention to the 
fact that a network analysis program is theoretically 
capable of converting a digital computer into a 
veritable analog computer of great versatility and 
power. While it may not be practical to write a 
single program to do this, at least two programs 
have been written for solving nonelectrical problems 
by network analytical methods: DYANA, 3 which 
handles either electrical or mechanical problems, 
and STRESS* which is based on a network approach 
to structural analysis. 


Computers are also influencing network theory 
by demanding methods of analysis adapted to the 
solution of computer-sized problems, As the 
engineer will appreciate, maximum "power trans- 
fer’' from the computer to his network problem 
requires a method which properly matches the 
computer to the problem. Traditional methods for 
“hand solution'' of networks are not necessarily 
best for use on a computer with networks of much 
greater size. The Laplace transform techniques 
fit this category and should at least be supplemented, 
if not supplanted, by numerical methods better 
adapted to the computer. One such method, 


described below, gives essentially all of the infor- 
mation obtainable by Laplace transforms. 


The objective of this paper is to highlight the 
importance of keeping the computer and its pecu- 
liarities in mind when developing the theory on 
which to base network analysis programs. Since 
matrix notation is the tool par excellence for des- 
cribing the network problem in a form suitable 
for programming, we shall use a matrix approach 
throughout the following discussion of various 
aspects of d-c, a-c, and transient analysis of net- 
works on a digital computer, 


Matrix Formulation of the Network Problem 


It is an interesting historical coincidence that 
just prior to the advent of the high-speed elec- 
tronic computer, the pioneering work of Gabriel 
Kron’ demonstrated convincingly the superior 
organizational powers of the matrix-tensor nota- 
tion in network theory. Kron emphasized not only 
the conceptual elegance of this notation but also its 
virtually autornatic way of handling the tedious 
bookkeeping chores of network analysis. He there- 
by laid some of the logical and procedural founda- 
tions necessary for programming a digital compu- 
ter to analyze networks — that is, both to compile 
and to solve the network equations automatically. 


The following matrix formulation of the net- 
work problem®»? is a direct outgrowth of Kron's 
work?:10 ang that of Roth?!-13 who provided an 
algebraic-topological explanation of the network 
problem as an abstraction. This formulation has 
been used at least in ae by such preerems as 
TAP, 14.15 yer -1,16 and ECAP,!? Kron's work 
has also been used as the basis for much of the 
computer analysis of networks encountered in the 
electrical power industry. 18 - 


From a physical point of view, the network 
problem is concerned with predicting the behavior 
of a system of interconnected elements in terms 
of the element characteristics and the manner in 
which these elements are interconnected, Viewed 
as a mathematical entity, the network problem 
involves the properties of an underlying topological 
structure, called a linear graph, and a superim- 
posed algebraic structure consisting of the inter- 
relations between certain quantities associated with 
the nodes, branches, and meshes of the graph. 


Obviously, the very act of interconnecting the 
network elements introduces certain constraints 
on the physical variables of the network problem 
and this gives rise to part of the bookkeeping chore. 
But the bookkeeping can be systematized in terms 
of matrices which describe the connectivity prop- 
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erties of a linear graph. Here is where topology 
is introduced, albeit in a very rudimentary form. 
Finally, the physical variables of the network 
problem correspond to the quantities which enter 
into the algebraic structure referred to above. 

The matrix approach nicely correlates these two 
points of view and has the additional merit of 
identifying precisely the separate roles of the 
element characteristics and of the interconnections 
in determining the overall behavior of the network. 


The ABCD's of Network Topology 


There are four topological matrices relating 
to a directed (or oriented) linear graph which are 
useful in nctwork theory. We shall describe them 
briefly. (More detailed descriptions are given in 
references 5 and 9), 
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Figure 1. Linear graph. 


Consider the graph of Figure la for which the 
branch-node connections may be tabulated as 
follows: 


Branch Initial Node Final Node 


1 
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This connection table, which is easily obtained 
from input data describing the network configura- 
tion, is a convenient starting point for compiling 
the network equations. (Usually, the nodes are 
numbered for the computer input, but letters are 
used here to distinguish nodes from branches. } 


An alternative way of storing this connectivity 
information in a computer is in matrix form as 
displayed in Figure lb, where the columns show 
which branches are connected to each node and the 
polarity of these connections. The way in which 


this A matrix is obtained from the connection table 
is obvious by inspection. (Branch 7, being short- 
circuited, requires both a +l anda -1 entry in 
column C so that both entries cancel.) 


Since each row of the A matrix contains both 
a +l and a -1l as its nonzero entries, its columns 
are linearly dependent — that is, their sum is zero 
so that any column is equal to the negative sum of 
all the other columns, Hence, any one column 
may be deleted since it contains redundant infor- 
mation. We choose to delete the column corre- 
sponding to node Ein Figure 1 and regard this 
node as the datum or ground node. The resulting 
A matrix, called the branch-node matrix, is 
shown in Figure 2b, 
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Figure 2. Tree-link structure of graph. 


If we choose a tree in the graph—for example, 
branches 1,2,3, and 4 shown as solid lines in 
Figure 2a—and classify the remaining branches 
as links, we may partition the A matrix into the 
submatrices Ap and Aj, referring to tree-branches 
and links, respectively. (An algorithm for per- 
forming a tree-link sort on the computer is given 
in Appendix Il of reference 15. The sequence of 
branches may need to be rearranged before parti- 
tioning the A matrix, but this is a simple task. ) 


It is easily shown that the submatrix Ap is 
square and has an inverse which we may determine 
directly from the graph of the tree in terms of the 
node-to-datum path matrix Bey for the tree. 9 


For the tree shown in Figure 3a, the columns 
of By, shown in Figure 3b, indicate which branches 
are included in the paths fram each node to the 
datum node. The desired relation between By and 
the inverse of Ap is: 


a1 sagt a) 


Ay T 


where pt is the transpose of the B,, matrix. (An 
algorithm for determining the By matrix from the 
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A matrix is given in Appendix Ill of reference 15.) 


By MATRIX 
NODE 
BRANCH A 8B i) 
1 +1 0 08 0 
2 o Oo -1 ie] 
3 7 1 1 0 
4 0 Oo 0 1 
(a) {b) 


Figure 3. Tree of graph. 


The next topological matrix of interest is the 
branch-mesh or circuit matrix, C, shown in 
Figure 4b. The columns of this matrix show which 
branches are included in each mesh. 
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Figure 4, Basic meshes of graph. 


If we adopt the convention that each link, together 
with its contiguous tree-branches defines a basic 
mesh whose orientation coincides with that of its 
defining link, then each basic mesh will contain 
only one link, oriented positively, as shown in 
Figure 4a, Consequently, the submatrix Cy;, 
referring to the links, is a unit matrix as shown 
in Figure 4b and it need not be stored explicitly in 
the computer. Only the submatrix Cr, which 
gives the path-in-tree for each basic mesh, is 
needed and this may be obtained as follows. 


It is a standard theorem that for any linear 
graph, the branch-node and branch-mesh matrices 
obey the fundamental relations, 


atc = 0 and/or ca = 0 (2) 


where the product is a null matrix. As a conse- 
quence, we may write the first of these expres- 
sions in partitioned form, 


t ot _ at t | 
[Aes] Cp] =A, Cp t AL 0 (3) 


We may solve this equation for C-»p by inverting 
A and using equation (1) to obtain the relation 


: t 
Ct Pa, i 


The final topological matrix we need is the 
cutset matrix D. Just as each link, together with 
a certain set of tree branches, defines a basic 
mesh, so each tree branch together with a certain 
set of links defines a basic cutset. This is shown 
in Figure 5a, 
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Figure 5. Basic cutsets of graph. 


The columns of the D matrix show which 
branches are included in each cutset. With each 
tree branch positively oriented in its correspond- 
ing basic cutset, the submatrix D7, referring to 
the tree branches, is a unit matrix while the sub- 
matrix Dy, referring to the links, turns out to be 
the negative transpose of Cy. (See reference 9, 
for example.) This is shown in Figure 5b, As a 
resuit, it follows that 


t t 
2 5 
ie - AL] BpZ AB, (5) 
L T er 


so that the cutset matrix D is readily computed in 
terms of A, By, and/or Cp. Except for pointing 
out that 
t t 

DC = Oand/or CD = 0 (6) 
this completes our description of topological mat- 
rices. Their use in performing the bookkeeping 
tasks in network analysis will be explained in the 
next section. 
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The Electrical Network Laws 


Consider as a typical branch in an electrical 
network the configuration shown in Figure 6, 


Figure 6. R-th branch of an 
electrical network, 


Here, the r-th branch is shown consisting of an 
impedance (or admittance) element Z,, (or Yp;) 
in series with a voltage source E,, and in parallel 
with a current source l,, either of which sources 
may have any value. Thus, there are three dis- 
tinct voltage and current variables to identify in 
this branch. As shown in the drawing, branch 
current i, enters the branch at its initial node but 
the element current J, actually passing through : 
Zr, is, according to the sign conventions depicted, 

J. e ys + *y (7) 
Similarly, although the branch voltage, measured 
from initial to final node, is e,, the element voltage 
V, actually appearing across 2, is 


V_ =E_ te (8) 
r r r 
with the sign conventions taken as shown. 


Clearly, since E, and I,, may be nonzero, Ohm's 
law for this branch must be written in terms of the 
element voltage and element current variables Vy 
and Jy. Accordingly, we may express Ohm's law 
for the entire network as the matrix ‘equations 


V=ZI (9) 


u 


and/or 


J=YV (10) 


where the vectors V and J consist of the entire set 
of element voltage and element current variables 
for the network and where the admittance matrix 
Y is the inverse of the impedance matrix Z. The 
diagonal terms of Z (or Y) are the self-impedances 


(or self-admittances) of each branch, whereas the 
off-diagonal terms are the transimpedances {or 
transadmittances) between pairs of branches. If 
corresponding off-diagonal terms are equal, they 
correspond to mutual impedances or admittances. 
Active devices, or so-called "dependent sources," 
may be represented by off-diagonal terms in the 
Zor Y matrix. 


In order to account for the constraints on the 
branch voltages e and branch currents i imposed 
by the interconnections, we invoke Kirchhoff's 
voltage and current laws. In their original form, 
these laws correspond to the matrix relations 


t 


Ce =0 (11) 


and 


(12) 


where the topological matrix fey acts as an operator 
which sums all the branch voltages around each 
basic mesh while the matrix A*t sums all the branch 
currents leaving each node. 


Although we may not always be aware of the 
fact, we are actually using alternative and mathe- 
matically equivalent forms of Kirchhoff's voltage 
and current laws when we say, ''the branch voltages 
are a linear combination of the node voltages" and 
"the branch currents are a linear combination of 
the mesh currents. '' These two statements corre- 
spond to the matrix equations, 


(13) 
and 
(14) 


where e' is the vector of node-to-datum voltages 
and i' is the vector of mesh currents. Here, as 
in equations (11) and (12), the topological matrices 
automatically do the bookkeeping for us. Inciden- 
tally, it is easy to derive equations (11) and (12) 
from equations (13) and (14) by premultiplying the 
latter by ct and at, respectively, and then invok- 
ing equation (2), The reverse derivations can also 
be made, proving that both forms of Kirchhoff's 
laws are equivalent. 


The basic interrelations between the network 
variables and the topological matrices may be 
summarized by a transformation diagram, given 
by Roth, 1 os shown in Figure 7. 
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(NODE) 


Figure 7. Algebraic structure of 
the network problem. 


This diagram, whose relation to equations (9) 
through (13) is evident, characterizes the alge- 
braic structure associated with a linear graph. 
The only additional relations are CtE=E, which 
defines the equivalent mesh voltage sources E', 
and Atl=I', which defines the equivalent nodal 
current sources I', 


By a proper interpretation of this algebraic 
structure, it can be shown that a wide variety of 
nonelectrical systems—e. g., mechanical and 
structural systems—fit into this algebraic pattern 
and so can bé treated as network problems. More- 
over, by extending this algebraic structure to 
characterize topological structures containing 
surface and volume elements as well as points 
and lines, a direct correspondence with the oper- 
ational structure of the vector calculus emerges, al 
Asa result, the use of network analogies for a 
wide class of partial differential equations — 
including Maxwell's equations—can be justified. 
The exploration and exploitation of this interesting 
topic is one of the broader aspects of network 
theory whose development the computer is helping 
to stimulate. 


The Electrical Network Problem 


A formal statement of the electrical network 
problem may be expressed as follows: "given a 
network whose configuration determines the topo~ 
logical matrices A and C, given its impedance 
matrix Z (or admittance matrixY), and given the 
voltage source vector E and current source vector 
I, find the branch voltages e and branch currents 
i such that equations (7) through (12) hold true. "' 


We shall discuss four distinct ways of attacking 
this problem. In each case, it will be necessary 
to introduce auxiliary variables such as the node 
voltages e' or mesh currents i' appearing in 
equations (13) and (14). In each case, it will also 


be necessary to invoke both of Kirchhoff's laws 
and, since each of these laws has two equivalent 
forms, the use of one form is tantamount to invok- 
ing the other. 


The Mesh Method. : 
dance form of Ohm's law given in equation (9), we 
introduce equations (7) and (8) and rearrange terms 
to obtain the expression 


Starting with the impe- 


(E-Z]) te = Zi (15) 
If we prernultiply by fou and take into account 
Kirchhoff's voltage law in the form of equation (11), 
we find that 
t te 

C (E-ZI) = C Zi (16) 
We may then incorporate Kirchhoff's current law 
by using equation (14) and writing the matrix equiv- 
alent of the mesh equations 


t t a 
C (E-ZI) = C ZCi (17) 


where ctzc is the mesh impedance matrix. 


If we solve this system of equations for the 
mesh currents i'—or write the formal solution for 
i' as 


i = (c'zc) “'ce-zn (18) 
where (ctzcy} is the mesh solution matrix--we 
may compute the branch currents i by using equa- 
tion (14) and finally obtain the branch voltages e 


by reverting to equation (15). 


The Node Method. If we start with the admit- 
tance form of Ohm's law given in equation (10), we 
may proceed to treat the network problem via the 
node method, following steps quite similar to those 
just described for the mesh method. Without des- 
cribing the steps in detail, we simply list the 
corresponding equations: i 


(I- YE) ti= Ye (19) 
A‘(1- YE) = a'ye (20) 
AL- YE) = A’yAe! (21) 
et = (atyay tat. YE) (22) 


In equation (21), A'YA is the nodal aera thelr 
matrix while in equation (22) (atya)-! is the nodal 
solution matrix, The desired variables e andi 
may be computed from equations (13) and (19). 


The Cutset Method. . In both of the foregoing 
methods of analysis, the auxiliary variables i' and 
e' constitute a "basis" for all the branch currents 
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ior branch voltages e, respectively. A basis 
such as e' is a linearly independent set of voltages 
in terms of which all the branch voltages may be 
computed, using equation (13). Any other basis 
besides the node-to-datum voltages will serve as 
well. In particular, we may use the set of tree- 
branch voltages e-; which are linearly independent 
and in terms of which the branch voltages may. be 
computed using the expression BS 


e=De 


" (23) 


Here Dis the cutset matrix. This is a statement 
of Kirchhoff's voltage law and it is the counter- 
part of equation (13), The counterpart of equation 
(12) for Kirchhoff's current law is 


t 
Di=0 (24) 


which follows directly from equations (6) and (14). 


Using these relations, the steps involved in the . 
cutset method of analysis are readily seen to be 


(I- YE)ti= Ye (25) 
p' (I- YE)=D'Ye (26) 
pir- YE)=D'YDe,, (27) 

_ (ptyp)")p‘U- YE) (28) 


Here, D, YD is the cutset admittance matrix and 
(DtyD)7! the cutset solution matrix. The variables: 
e and i may be evaluated using equations. (23) and 
(25). 


The Mixed Method. In addition to these three 
standard methods of analysis, there is a fourth 
which is beginning to be more widely used, namely 
the mixed method. This method, which is the 
basis for the writer's extension 5 of Bashkow's 
A-matrix formulation of the transient problem 
and for a new method of a-c analysis, © can also 
be used for d-c analysis. Moreover, it underlies — 
both Bryant's treatment” of Bashkow's sees 
tion and the so-called state variable approach 
to transient analysis. : 


The idea behind the mixed method is that Ohm's 
law may be expressed in a mixed form involving 
the admittance matrix for certain branches of the 
network and the impedance matrix for the rest of 
the branches. Thus, we may write Ohm's law in 


the form 
I ag i) Vv (29) 
y y y 
vil “jo 2 5 
PA Zz Z 


where the subscripts y and z refer to the branches 
which are treated as admittances or as impe- 
dances, respectively. With this equation, or its 


equivalent, 
L. $i Y 0 E te (30) 
y a y ¥ ¥ 
E +e 0 Zz tt 
zZ zZ zZ Zz Zz 


as a starting point, we may analyze the admit- 
tance branches by the cutset method and the 
impedance branches by the mesh method--virtu- 
ally considering the two sets of branches as 
comprising separate subnetworks. We then 
couple the resulting two sets of equations together 
by means of an appropriate topological matrix. 
The procedure is as follows. 


Taking all the y-branches (admittances) first, 
a tree-link sort is performed. Any y-links thus 
identified will define basic meshes which contain 
only y-tree branches since no z-branches (impe- 
dances) have been considered up to this point. 
We will designate the paths-in-tree for these y- 
meshes by the submatrix CTyy. 


The tree-link sort is then resumed, taking the 
z-branches into account. Now, however, the z- 
links define basic meshes whose paths-in-tree 
may involve both y-tree branches and z-tree 
branches, designated by the submatrices Cryz 
and Cyrzz, respectively. Thus, the Cy matrix 
for the entire network has the form 


(31) 


If we express the y-branch voltages e, in 
terms of the y-tree branch voltages, eTy> and 
recall equations (5) and (23) we have 


ae as = 2 
y ery = Vey ery Beaty (32) 
t 
e -c 
Ly Tyy 
where D, is the cutset matrix for the y-branches 


alone with the z-branches regarddd as being open- 
circuited. For the entire network, the corre- 
sponding relation is 
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Uny 1?) ony 
g Ure has 
= t a3 
Sayy a 
ec Tyz . Cree 


where the D matrix for all of the y- and z-branches 
is required. 


Similarly, we may express the z-branch 
currents ig as a linear combination of the z-link 
currents ine which are identical with the mesh 
currents, obtaining the relation 


(34) 


Here C, is the circuit matrix for the z-branches 
alone with the y-branches regarded as being 
short-circuited. Again, for the entire network, 
we have 


Cc 
Ty Tyy 
i 0 
Tz 2 (35) 
ate Vy 0 
*Lz 7g Mie 


where the C matrix for all of the y- and z-branches 
is required, 


Now it is true for the entire network that 
Dti=0Qandthat Cte =0, Asa consequence of 
these relations and equations (32) through (35) it 
may be shown that 


t. fs 3 
= « = 36 
Dy y [Pry Cryy| “Ty CryatLe (36) 
and iby 
t fat - oct 
Cc e. om v,, | en, Cry 2’ Ty (37) 
*L 


Returning to equation (30), we may rearrange 
terms and write 


(l_-Y_E ¥y 0 e i 
y oy ylefy y y 


(E.-Z1)} Jo z i e 
zZ ZZ 


(38) 


The next step is to premultiply the first row of 
this equation by Dt and the second row by ct and 
to replace e 


y and i, by the expressions given in 
equations (32) 


and (34). The result thus obtained, 
t t 
Di -Y Ej] f[p'y po e Di 
yoyoy y yyy Ty| _ YY] (39) 
ties. t t 
C (FE -2ai Z i 
us A 2 2 0 So oe ha ee 


is —-except for the second term on the right — 
exactly what would follow from using the cutset 
method of analysis on the y~branches, with the 
a-branches open~circuited, and the mesh method 
on the z-branches, with the y-branches short- 
circuited. Hence the name mixed method. 


The final step is to replace the terms involving 
i,, and e, by taking advantage of the relations de- 
veloped in equations (36) and (37). When this is 
done, we find that 


t t 
D (lL - YE ) DYD -C e 
yy yy. yy y Tyz TY] (40) 
t t t . 
oe a) Coys cr ar “Lz 
where the topological matrix C effectively 


5 z 
couples the two systems of cutset and mesh equa- 
tions. 


Methods of Solving the Network Problem 


The foregoing discussion relates to four 
general methods of network analysis -- that is, 
"separation of (the network problem) into con- 
stituent parts or elements" (Webster's Collegiate 
Dictionary). The present section refers to methods, 
primarily numerical, for solving network problems 
of three main types: d-c, a-c, and transient, 


D-C Network Problems 


The most straightforward way of solving d-c 
network problems is by Gaussian elimination2> of 
the variables in the mesh or nodal equations. The 
amount of computation required is about n3/2 
operations (multiplications and additions) where n 
is the number of variables. If several solutions 
are desired for a given network with a succession 
of voltage and/or current sources, the procedure 
generally followed is to compute the mesh solution 
matrix or nodal solution matrix, at a cost of about 
n3 operations, and then to multiply each different 
source vector by this matrix, at a cost of n®@ oper- 
ations per vector. 


The matrices involved in most network prob- 
lems are so sparse, however, that it is worthwhile 
to take advantage of this sparsity in programming 


the actual method of solution. Unfortunately the 
electronic practitioners of the network analysis 
programming art have not been as alert to this 
possibility as their confreres in the electrical 
power industry. This may be due to the fact that 
the power engineer is being forced to consider much 
larger problems~-in the order of 1000 to 2000 
variables--than is the electronic engineer, for 
whom 50 to 100 variables seems to be an upper 
limit at present, 


In a significant paper on this topic, Sato and 
Tinney26 have described two techniques for taking 
advantage of the sparsity of the nodal admittance 
matrix whereby the number of computations can be 
greatly reduced at the same time that the require- 
ments for computer memory are decreased. These 
authors show that the order in which the variables 
are processed in the usual Gaussian elimination 
scheme strongly affects the number of computations 
involved in treating a sparse matrix. This question 
has also been considered by Parter. at 


One rule which considerably reduces the com- 
putation is to eliminate, at each stage, that variable 
for which the corresponding row in the still-to~be- 
triangularized matrix has fewest nonzero elements, 
Although this simple rule cannot produce an optimum 
strategy, it does yield improvements in the order 
of 4 to 1 in computing speed. 26 This is because 
the rule helps to preserve sparsity in the triangular 
matrix which results from the elimination procedure. 
A better rule, proposed by Bree, 28 is based on a 
graph theoretical representation of the problem 
similar to that described by Parter. 


Sato and Tinney describe a second technique 
which not only reduces storage requirements, by 
eliminating the need to store the inverse matrix, 
but also produces a codified structure that can be 
used in place of the inverse matrix. This structure, 
which preserves the elimination operations, can be 
used to transform a "source vector" into a "solution 
yector'' much more efficiently than can be dane by 
matrix-vector multiplication using the explicit 
inverse matrix, which is hardly ever sparse. 


As Sato and Tinney point out, the Gaussian 
elimination (triangularization) procedure and the 
back-substitution process can be represented as a 
product of elementary matrices. These are mat- 
rices which differ from a unit matrix by either one 
diagonal or one off-diagonal entry. The essen- 
tial entry of each of the elementary matrices for 
the Gauss elimination process can be stored on or 
below the diagonal of the resulting triangular matrix; 
the essential-entries for the back-substitution pro- 
cess are simply the superdiagonal terms of this 
matrix, Thus, the space occupied by the original 
matrix can be overlaid by this codified structure 
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which is equivalent to a "product form of the 
inverse, "29 For very sparse matrices, special 
storage schemes for compacting this information26 
result in great savings in both memory space and 
computation time. Even for a full matrix, com- 
puting the product form of the inverse requires 
about half the operations needed to compute the 
inverse explicitly. 2 


The order of elimination of the variables 
affects not only the amount of. computation {in a 
sparse matrix) but also the accuracy of the com- 
putation--in any matrix. To minimize round-off 
errors, the rule generally followed is to inter- 
change rows and columns, at each stage of the 
elimination procedure, so that the largest possible 
term appears in the "pivotal" position. 25,29 (This 
is the position on the diagonal corresponding to the 
next variable to be eliminated.) In many network 
problems, the matrices are diagonally dominant, 
that is the diagonal term exceeds the sum of the 
magnitudes of the off-diagonal terms in the same 
row or column. Even so, for best accuracy, the 
interchanging technique is worthwhile and not too 
expensive in machine time. A combination of this 
technique with that of Sato and Tinney should be 
feasible. 


Sometimes, inaccuracies arise which are not 
related to the round-off errors introduced in the 
process of solving equations or inverting 2 matrix. 
For example, suppose the node method is being 
used and it is desired to compute the current in a 
branch whose node voltages are almost equal. 
Clearly, a severe loss of accuracy will occur when 
these two node voltages are subtracted during the 
evaluation of equation (13), even if e' is known to 
full machine precision. This kind of problem, 
then, is a function of the source vectors E and L 


One way of avoiding this difficulty, which was 
actually used in TaP,! is to compute currents 
using the mesh method and voltages using the node 
method. An alternative which has not been tried, 
but merits investigation, is the mixed method. 


Kron's method of interconnecting solutions 
provides some tools which are very useful for 
certain applications. In particular, if it is desired 
to study the effect of varying a single resistance, 
letting it take on a succession of values, Kron's 
method provides a much more efficient way of 
handling this problem than solving each case from 
scratch, This method, which is the basis of the 
MODIFY feature of ECAP, permits the nodal solu- 
tion matrix for the original network to be updated 
in accordance with changes in one or more branch 
resistances. The amount of computation involved 
is trivial compared to that of solving the nodal 
equations or inverting the nodal admittance matrix 
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8,9,10 


all over again. 


A simple way of understanding the use of 
Kron's method in this application is the following: 
Consider the network configuration shown in 
Figure 8a which consists of a single link connected 
from node p to node q of a tree all of whose branches 
converge on the datum node. The path-in-tree 
for the single mesh in this network is given by the 
Cy matrix denicted in Figure 8b. 


0 
i?) 
° 
-1 |p-th row 
Cre ° : 
o 
1 |q-th row 
° 
0 
(a) (b) 


Figure 8. Connection of single link to a tree. 


. P t % 
The mesh impedance matrix C ZC for this net- 
work is easily shown to be 

t 


c'zc=cl,z 


peop tt 


L (41) 
where Zp is the impedance matrix for the tree 


while Zy, is the impedance "matrix'' of the single 
link being added to the tree. 


Now by substituting equations (14) and (18) into 
equation (15), it can be shown that 


e=(z- zc (ctzc) tctz] (1-YE) (42) 


Then, by partitioning this equation into submatrices 
and subvectors relating to trees and links, it can 
further be shown that? 


: to -1 at t 
e,=[Z CC ZC) CZ, ]D(I-YE) (43) 


pleg 2p 
Finally, by comparing this result with equation 
(28) and using equation (41), it follows that 
t -1_t 


Lp CglCpZpCpt Zp) CpZep (44) 


t -1 
(D'YD) =Zn 
Thus, we have derived from the mesh method 

an expression for updating the cutset solution 
matrix for any tree to which a single link has been 
added. (Remember that for a tree, the impedance 
matrix Zp is identical with the tree solution matrix 
qptyD)-} since D=D,, is a unit matrix.) But in 


view of equations (4) and (5) we can extend this to 
an expression for updating the nodal solution ma- 


trix as well, ? namely 
t -1 t 308 -1 
eZ eZ A 
{A YA) a Zs AZ 4, t4,) AL, (45) 


t F 

where Z)*B,.2.B i 

Returning now to our original problem, that 
of updating the nodal solution matrix when a single 
branch resistance (between nodes p and q) is. 
changed, we assume that this change is effected 
by adding another resistance in parallel, either 
positive or negative in value, We assume also 
that the nodal solution matrix for the original net- 
work has been computed, This matrix may then 
be interpreted as the impedance matrix Zq of a 
tree having the ''star'' configuration shown in 
Figure 8a. Since Br for this tree is a unit matrix, 
it follows that we may substitute Zp for Z, in 
equation (45) which thus enables us to update the 
nodal solution matrix one link at a time. 


Equation (45) may be used recursivel for 
adding several links, one after the other.” If this 
technique is used, however, to compute a sequence 


of solutions as a given resistance is stepped through 


several different values, round-off errors will 
accumulate from one step to the next, If the range 
of variation of resistance is large enough and many 
steps are taken, the errors may eventually become 
serious enough to vitiate the results. A better 
technique, but one requiring more memory space, 
would be to update only the original nodal solution 
matrix each time. This would preclude the accum- 
ulation of errors from one step to the next. 


By suitably extending this idea of updating the 
nodal solution matrix to take into account changes 
of a single resistance yalue, differential changes 
in resistance can be dealt with. In other words, 
the partial derivative or sensitivity of the nodal 
solution matrix with respect to any resistance 
change can be computed explicitly. 


Rather than describing how this may be done 
starting from equation (45), we may derive a 
similar and broader relation by differentiating the 
nodal equations themselves. Thus, from equation 
(21) it follows that 


Atidl-¥ db-dY E) = A‘ay Ae! +A’ YAde! (46) 
or 
de' =(Atyay”! a! fal-¥ dE-d¥ (Ete)] (47) 


The latter expression gives the variation in 
node voltages as a function not only of changes in 
resistance (specified as d¥) but also of changes in 


the Land E vectors. This equation is used in ECAP 
for computing sensitivities. 


The question of solving nonlinear d-c network 
problems is an important one and difficult to answer 
with complete assurance. One standard technique 
that is useful in many cases, however, is the 
Newton-Raphson iteration method. 30 For solving 
a single nonlinear equation f(x)=0, this method uses 
the recursion formula, 


1-H)! £05) (48) 


itl 
where f* (x;) is the value of the derivative df/dx at 
the point x 


This formula may be extended to handle a 
system of n nonlinear equations, 
(49) 


het Eee Parra x,)=0; jal,2,...,n 


2 
by generalizing the role of the derivative. If we 
define the Jacobian matrix H in terms of its jk-th 
element, 


Hy = 06/8 x, (50) 
then equation (48) may be replaced by 
C -1 
*4 5%, -H f. (51) 


where x; and f. are vectors, An adaptation of this 
method has been used with some success in Tapi5 
and in the NET-1 program, !6 


When the Netwon-Raphson method converges, 
it does so quadratically. Thatis, at each iteration 
the error is proportional to the square of the pre- 
vious error. However, the method may diverge 
under certain circumstances, for example, in the 
case of transistor switching circuits with positive 
feedback. 


In a recent paper, Sroyden?= has described a 
class of methods which promise to remedy two of 
the main drawbacks of the Newton~-Raphson method, 
namely: the need to evaluate and invert the Jacobian 
matrix and the liability to diverge. To the writer's 
knowledge, Broyden's method has not yet been 
applied to nonlinear network problems, but it de- 
serves to be investigated. 


Another approach to solving the nonlinear net- 
work problem has been described by Katzenelson. 
This method converges in a finite number of steps 
if the resistance characteristics are continuous, 
monotonically increasing and piecewise linear. 
Under these restrictions, the network problem has 
a unique solution, The solution is found by solving 
a succession of linear network problems which are 
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locally equivalent to the actual problem. 


This sequence of linear problems is generated 
by gradually increasing the voltage and/or current 
sources from zero to full value and solving each 
locally linear problem within 2 finite region whose 
boundaries are determined by the breakpoints in 


2 


the piecewise linear characteristics of each resis- 
tor. Thus, a solution curve, consisting of a finite 
number of straight line segments, is traced out in 
the multidimensional space of solutions, ending at 


the desired solution point. 


The key to the computational effectiveness of 
Katzenelson's algorithm is the use of Kron's 
method for updating the cutset solution matrix by 
means of the link-at-a-time technique indicated in 
equation (44). Although the method might be used 
successfully for nonlinear problems with multiple 
solutions (e.g, tunnel diode circuits, multivibrators, 
etc.) it is guaranteed to converge only under the 
conditions mentioned which insure a unique solution. 


A-C Network Probleins 


Since all of the theoretical developments up 
through cquation (47) apply to complex as well as 
to realnumbers, steady state a-c network problems 
may be handled using any of the four basic methods 
of analysis described above as well as Kron's 
method, For example, by programming the Gaus- 
sian elimination procedure for complex arithmetic, 
the mesh or nodal equations for any a-c problem 
may be solved directly, at a computational cost of 
approximately four times that for real arithmetic. 
However, since the coefficient matrix is frequency- 
dependent, the equations must be solved or the ma- 
trix inverted afresh at each new frequency. Aside 
from the computational cost, this approach is 
usually quite reliable, although it may involve 
reduced numerical accuracy in the vicinity of a 
resonant frequency. 


Anew approach ta the a-c problem, 6 which is 
based on the mixed method of analysis, tackles the 
job of matrix inversion by first solving the much 
more difficult problem of computing the eigen- 
values and eigenvectors? of the matrix, At first 
glance, solving the eigenvalue problem appears to 
be "doing it the hard way," since the amount of 
computational effort required is an order of mag- 
nitude greater than that for inverting a matrix. 
This would indeed be the case if only a single in- 
verse roatrix were needed. But in calculating 
frequency response, a new inverse, or more 
precisely, a new solution vector, is needed at 
each frequency. In this case, using the eigen- 
values and eigenvectors turns out to be just the 
right thing to do because it reduces the problem 
to one of inverting a diagonal matrix at each 


new frequency. 


Starting with equation (40) of the mixed method, 
we can derive the differential equations for an RLC 
network in first order form as follows: All capaci- 
tors are treated as admittance branches and hence 
may have a nonzero conductance in parallel. Thus 
the admittance matrix becomes 

Y =pK+G (52) 

¥ 

where p=d/dt and the symbol K is used for capaci- 
tance so as not to conflict with the symbol C for 
the circuit matrix. On the other hand, all inductors 
are regarded as impedance branches with (possibly 
nonzero) resistance in series so that the impedance 
matrix is 

Z =pLtiR (53) 
{Branches consisting of resistance only can also be 
handled. But since these give rise to algebraic 
equations!5 which unduly complicate the discussion, 
we shall not consider them here, ) 


Using equations (52) and (53) together with 
equation (40), we may write the differential equa- 
tions (for constant K and Li) in the form 


D'cD 


ct 


Tyz 


DKD 0 -c 


Tyz 
crc 


pe 
TWh (54) 


t ‘ : 
(0) Cc Lc pip i L 
where the symbols Ip and E, designate the driving 


currents and voltages. 


In more compact form, equation (54) may be 

written as 

Px + Qx = f(t) (55) 
where x = dx/dt and x is the state-variable vector 
of capacitive tree-branch voltages and inductive 
link currents. By inverting P we may rewrite 
equation (55) in the form 

x = Ax + g{t) (56) 

7 -l 

where g=P le and A=-P Q, which corresponds 
to Bashkow's A matrix. %* (This matrix is not to 


be confused with the branch-node matrix A, which 
will not be referred to any further. ) 


For steady state a-c problems, the driving 
function may be written as 


g(t) = ae (57) 
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with s a complex number, 
sponse then is 


The steady state re- 


st 
x(t) = xe 


(58) 
substituting these relations in equation (56) we 
have the linear system of equations 


(sU-A)x, = Bo (59) 


where g_ is a known vector and x_ is to be deter- 
: ° 
mined, 


The formal solution for equation (59) is 


-1 
Pa -A 
x, = (sU-A) gg, (60) 
and so we need to find an effective way af compu- 
ting (sU-A)~ 1, This may be done using the cigen- 
values A; and eigenvectors yj; of the A matrix as 
follows. 


Let Abe the diagonal matrix of eigenvalues 
and X the matrix whose columns are the corre- 
sponding eigenvectors. The eigenproblem, then, 
is characterized by the equation 


AX = XA (61) 
which may be recast into the form 
x lax = A (62) 


This defines a similarity transformation”? which 


diagonalizes the matrix A, 


The same transformation also diagonalizes 
(sU-A) giving 


x7 (su-a) X= sU-A (63) 


Solving this expression for (sU-A) and then invert- 
ing, we find that 


1 


(au-a)7) = x (st-Ay xX” (64) 
so that equation (60) becomes 
-ll-1 
= - 65 
x, = X(sU-A) X gy (65) 


This is the basic result of the new method, requir- 
ing the inversion of the diagonal matrix (sU- A) 
at each frequency. 


E we letd = x g_, which we may calculate 
once for all, then the i2th element of the Xy vector 
may be computed using the relation 


no X..d. 
nie SS) (66) 
oe APP ee A 


This partial fraction expansion of the state variable 
x, May be evaluated at a computational cost pro- 
portional to n operations per frequency as opposed 
to n? operations if equation (59) were solved by 
Gaussian elimination, 


Now the cost of solving the eigenproblem is 
not negligible since it is equivalent to solving equa- 
tion (59) at roughly 10 frequencies. Accordingly, 
the present method would not pay off if solutions 
were required at only a few frequencies, But for 
a frequency response calculation, this method is 
very effective. 


It is quite easy to extend this method to the 
computation of driving point and transfer functions 
of the network and to derive sensitivity formulas, 
For example, the frequency sensitivity of x is 


1 


ee isU AL £, (67) 


(See reference 6 for details.) But one of the more 
interesting features of the method is the fact that it 
is far better suited for computer-size problems 
than the traditional Laplace transform techniques 
involving ratios of polynomials and the poles and 
zeroes thereof. 


In particular, the task of computing the coeffi- 
cients of the polynomials in a network function 
P (s) / Q (s) is not only time-consuming but also 
prone to serious numerical inaccuracies, especially 
when the polynomials are of high degree. The so- 
called “topological formula" approach to computing 
these network functions involves finding all the 
trees of a network and then computing the sum of 
the corresponding tree-admittance products. But 
the number of trees may run into millions for a 
network with only 20 nodes and 40 branches, 2 And. 
even if this were not enough of an impediment, the 
computation of the roots of the polynomials P (s) 
and Q {s) is hazardous because these roots may be 
extremely sensitive to errors in the coefficients. 25 


In the writer's judgment, therefore, the poly- 
nomial approach is just not matched to the network 
analysis tasks which the computer is called upon 
to handle. The eigenvalue approach is much better 
suited and gives all of the theoretical information 
that the Laplace transform methods are designed 
to provide. For example, the eigenvalues are 
identical with the poles of the network functions, 
as equation (66) indicates. Moreover, any network 
function desired may be computed straightforwardly 
and its sensitivity obtained, either with respect to 
frequency or with respect to any network parameter. 
Finally, even the pole sensitivities can be calculated. © 
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This is not to say the eigenvalue method is 
free of all difficulties. Even with the Q-R method 
of Francis34 which, according to Wilkinson, 25 is 
"the most effective of known methods for the solu- 
tion of the general algebraic eigenvalue problem," 
it is possible that numerical errors of a crippling 
nature may arise when networks having a large 
ratio of greatest to smallest eigenvalue are treated. 
The difficulty with this kind of "pathological" net- 
work is that although the large eigenvalues can be 
computed accurately, the small ones cannot. But 
it is the small ones which are principally respon- 
sible for determining the behavior of the network, 
as equation (66) implies. 


One possible remedy for this situation is to 
solve the eigenvalue problem for the inverse 
matrix Av}, Then the situation is reversed, for 
the large eigenvalues of a7}, which will be obtained 
accurately, are the reciprocals of the small eigen- 
values of A. Computing acl may be error-prone 
in its own right because of the large ratio of eigen- 
values. But if the matrix Q in equation (55) is 
nonsingular, then we may use the relation 


-1 1 


A" =-Q°°P (68) 


to compute a more accurate inverse for A. This 
question remains to be investigated. 


Fortunately, the problem of multiple eigen- 
values presents no difficulty since the Q-R method 
is able to handle this situation quite effectively. 

If, however, the matrix is such that it does not have 
a complete set of linearly independent eigenvectors, 
it may be necessary to perturb the matrix slightly 
so as to avoid this difficulty, From an engineering 
point of view, this subterfuge would probably be 
completely acceptable in most instances, 


Transient Network Problems 


Of the three types of network problem under 
discussion, transient problems present the greatest 
computational difficulty. This is because the 
numerical integration of the differential equations 


is burdensome for nonlinear networks and extremely 


slow when the network has a wide spread of eigen- 
values. The nature of these problems will be 
discussed presently. 


Analytic Solution of Linear Problems. To 


fix ideas, let us first consider the analytic solution 
of the state-variable equations for a linear network 
with constant parameters. Referring to equation 
(56), we assume a solution of the form 


x(t) =e fy(t) (69) 


where v(t) is a vector function of time and where 


At, x : F 
e is defined as the matrix function 


(ag? , (as? 


Al 
e floaty AL 
! 3! 


tone (70) 


By differentiating equation (69) and using the defi- 
nition of e“t it follows that 


At 


a f At- At. 
x=Ae t t 


vte v=Axte v (71) 


Combining this with equation (56), we obtain a 
differential equation for v(t) in the form 


(72) 


a 


See 


whose solution may be written as 


h 
v(t th) = v(t) + ie Peal ey +7 dt (73) 


4 
Multiplying this equation by et Bl aa using equa~ 
tion (69) we find that the solution to equation (56) is 
Ah h A(h- 
x(t th) =e" x(t) f a Tht +7)dT (74) 
° 


This solution is exact, but its utility depends on the 

computation of eAh, 

.. 2,34 suns 
Liou has proposed the use of a finite number 

of terms in the series expansion of e according 

to equation (70), asserting that this method ''com- 

pletely avoids eigenvalue problems. "' In one sense, 

this assertion is true in that Liou's method does 

not require the calculation of eigenvalues. On the 

other hand, the method runs ‘headlong into another 

kind of eigenvalue problem which seriously limits 

its usefulness; namely, that when the matrix A has 

a large eigenvalue, the integration step size h must 

be kept small, roughly less than the reciprocal of 

|AmaxtA)| , in order to permit rapid convergence 

of equation (70). Precisely the same limitation 

arises in numerical integration using, say, the 

Runge-Kutta or Adam's method, where h must be 

kept small to avoid numerical instability. 30. More 

about this later. 


If, instead of avoiding the computational eigen- 
value problem, we solve it, then equation (74) be- 
comes as practical to use as equations (65) and (66). 
From equation (70) it can be shown that the trans- 
fromation which diagonalizes A also diagonalizes 
eAh, yielding the relation 


xl Aby . eft (75) 


h, - ‘ ih 
where eM is a diagonal matrix whose terms ent 
are trivial to compute. Solving this expression for 
eh, we may recast equation (74) into the form 
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Ah 


x(tth)= Xe xe Mth-7) 


-i 
X x(t) + x” Taletride (76) 


and our computational problem reduces to one of 
evaluating the integral expression. If the driving 
function g (t+r) is approximated as a polynomial in 
7, the integration can be carried out exactly, term 
by term. This provides a very satisfactory way 
of evaluating the integral on a digital computer. 


Equation (76) is instructive from a theoretical 
For example, if g (t) = 0 and we 
apply initial conditions to the network which corre- 


point of view. 


spond to, say, the i-th eigenvector x y» then we 
may write the solution as 


x (t) = Xe At xly (77) 


i 
: -ll. , ; or re 
But since X Xis a unit matrix, X “yj is equal to 
the i-th column of a unit matrix. (Remember that 
the columns of X are tl the eigenvectors x;.} Thus, 
the eigenfunction eri' is singled out and equation 
(77) reduces to 
e Ait 


x (t) = (78) 


In other words, the i-th eigenvector prescribes 
the particular initial conditions required to excite 
the i-th characteristic mode (eigenmode) of re- 
sponse and suppress all the others. For arbitrary 
initial conditions, given by a vector x9, all of the 
eigenmodes will be excited, in general. But the 
matrix-vector product xo lx automatically pro- 
duces the appropriate linear combination of eigen- 
functions in terms of which to express the solution, 


For the case where g (t) # 0 the "initial condi- 
tions' are continually changing, so to speak, and 
so an additional term involving a convolution inte- 
gral of the eigenfunctions and driving functions is 
required, as in equation (76), 


Numerical Integration. The numerical methods 


of integration currently in use are, for the most 
part, based on the use of finite difference approxi- 
mations to the actual differential equations to be 
solved. For a single differential equation of the 
form 

dx/dt = f (x, t) (79) 
a solution may be effected by the general finite 
difference expression 


M M 
2S To Aes pha o£ bt. (80) 
n z 1 nel jn-j 

isl jz0 


where h is the integration interval, x,7* (nh) and 


fn-j= (yj. (n-j)h). The number and values of the 
coefficients a; and b; are determined by the number 
of terms of the Taylor expansion of x (nh) which are 


to be matched by the integration formula so defined. 


If bo=0, the integration formula is called a 
“predictor, since the new value of x, is predicted 
in terms of previous values (x,_;) and previous 
derivatives (fn-;) only. If b, #0, the formula is 
called a "corrector," since it is generally used in 
conjunction with a predictor formula to provide a 
corrected value of x,. This is accomplished by 
using the predicted value of x, to compute an approx~ 
imate value of the 'new" derivative f,, which is 
then used in equation (80). 


Now, in the simple case of one homogeneous 
linear differential equation with constant coefficient, 
where f (x, t) = g x (t), equation (20) may be written 
in the form 


(1-hgb.) x (a, thgb)) a ie (a, thgb x70 


(81) 


Here at least one of the coefficients a, or b, is non- 
zero and all coefficients of higher index are zero. 
That is, r=max (M,N). If we assume a particular 
solution of the form x,= p™, then equation (81).may 
be reduced to the polynomial equation, 


r r-1 , 
(1-hgb,)p - (a, thgb,)p meee ee (a, thgb_) iy (82) 


Obviously, any value of p which is a root of this 
characteristic equation is a candidate for the par~ 
ticular solution x,= p™, but the general solution 
takes the form of a linear combination of particular 
solutions, namely — 


(83) 


where the py, are all the roots of equation (82). 


One of these roots, lying near the point 1 + j0 
in the complex plane, will be the principal root 
since it primarily determines the desired approxi- 
mation to the solution of the differential equation in 
question. All the other roots are parasitic: roots 
and it is clear from cquation (83) that if any parasitic 
root lies outside the unit circle, x, will eventually 
diverge. This situation is called "numerical insta~ 
bility'' because it vitiates the solution. {Actually, 
Xy may diverge because the principal root is outside 
the unit circle, but this corresponds to a divergent 
solution to the differential equation itself. } 


In order to prevent numerical instability, the 
integration step h must be kept small enough ‘to 
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insure that all of the parasitic roots of equation 
(82) remain inside the unit circle. In the case of 
a system of linear differential equations with con- 
stant coefficients, as typified by equation (56), the 
requisite condition on h for most integration for- 
mulas is 


ata a t«< (84) 


Here \A\max (A) is the largest eigenvalue of the 
coefficient matrix A while the constant c is in the 
order of 1 or 2, depending on the integration 
method used. Even in the case of nonlinear dif- 
ferential equations, an equivalent restriction, 
involving the (instantaneously) largest eigenvalue 
of the Jacobian matrix, must be imposed to pre- 
vent numerical instability. 


It is just this constraint on the integration 
interval which makes most integration techniques 
so painfully slow in the case of network problems 
having a wide spread of eigenvalues. For the 
largest eigenvalue (or, equivalently, its reciprocal, 
the smallest time constant) controls the permis- 
sible size of the integration step, h. But the 
smallest eigenvalues (largest time constants) 
control the network response and so determine the 
total length of time over which the integration 
must be carried out to characterize the response, 
In the case of a network with a 1000 to | ratio of 
largest to smallest eigenvalue, for instance, it 
might be necessary to take in the order of 1000 
times as many integration steps with, say, the 
Runge-Kutta method as with some method which 
was free of the restriction given by equation (84). 
This so-called ‘minimum time-constant'' barrier 
is without doubt the biggest obstacle to be over- 
come in the transient analysis of networks. 


There are, however, some integration for- 
mulas in common use for which the constraint on 
h given in equation (84) does not apply; for example, 


x =X thf (85) 
n n-1 n 
and 
= + f +f 
ar ea Bi n at! 186) 


2 


Both of these have characteristic equations of the 
first degree and so have no parasitic roots. But 
the main disadvantage in using these formulas is 
that they are of such low order accuracy, particu- 
larly equation (85), that the integration interval 
must often be shortened anyway to insure sufficient 
accuracy. Another disadvantage is that both are 
implicit formulas because they involve f), which 


to equation (56), therefore, these formulas require 
the solution of a linear system of equations (or 
matrix inversion) which complicates their use. 


Schemes for retaining the numerical stability 
inherent in equations (85) and (86) without the need 
to solve simultaneous equations have been reported 
by Stineman3# and by Richards, Lanning, and 
Torrey, 35 but the resulting methods are again of 
low order accuracy. Treanor, on the other hand, 
describes a method of high (fourth order) accuracy 
which is identical with the Runge-Kutta method for 
short integration intervals but retains numerical 
stability even for very large integration intervals. 
His method, however, involves computation of the 
matrix function e“} and so presents no real advan- 
tage over the eigenvalue approach for linear prob- 
lems, In the case of nonlinear problems, the 
difficulty in using Treanor's method is compounded 
by the fact that e“h must be recomputed at each 
time step. 


Quasi-Analytic Solutions. A quite different 
attack on this problem has been reported by Cohen 


and Flatt37»38 and by Certaine. 3 (Incidentally, 
the title of Certaine's paper is misleading since 
what he calls "time constant" is the reciprocal of 
the time constant usually referred to by electrical 
engineers.) The essence of this attack is the 
following: in equation (56) the diagonal terms of 
the coefficient matrix A are handled analytically, 
as was A itself in equation (74), while the off- 
diagonal terms of A are lumped together with g (t) 
and considered as part of the driving function. 


In other words, letting A=D+B, where D con- 
tains all of the diagonal terms of A while B contains 
all the off-diagonal terms, we may write equation 
(56) as 

x=Dx + (Bx +g (t)] (87) 
with Bx + g (t) serving as a driving function. The 
analytic solution to this equation can be obtained 
by modifying equation (74) to read 


Dh h D(h-7t)- 
x(tth)=e~ x(t) +f, gotham) [Bx (ttr) + g(ttrJdr 
(88) 
Dh. : : é 1 
Here, ¢ is easy to compute, since D is a diagona 
matrix. The integration can be carried out by 


using a polynomial approximation of any desired 
degree for Bxtg. (The method of integration used 


in NET -1 is based on a modification of this approach.) 


When the matrix A is dominantly diagonal, 
the use of equation (88) may permit a much larger 
integration step size than that allowed in the usual 


cannot be computed until x, is known. When applied methods. However, for network problems, this 
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property of diagonal dominance of the A matrix 
does not generally obtain and so equation (88) is 
not of itself sufficient to break through the inte~ 
gration barrier. 


The writer has found that equation (88) has its 
own peculiar numerical instability which seems to 
be avoidable by choosing the integration interval 
subject to a restriction of the form 


plana (Bl sc (89) 


where cis a small constant. In other words, the 
largest eigenvalue of the off-diagonal matrix B 
controls the stability of the numerical integration 
process used in applying equation (88). 


From a theoretical point of view, it is possible 
to choose a similarity transformation on A of the 
form 

=Y a ee 

TAT =A=D+B (90) 
where Dis the diagonal component and B the off- 
diagonal component of A, such that Amax (B) is 
arbitrarily small. Indeed, if T is so chosen that 

x (B) = 0, then DD contains all the eigenvalues 
of A, (In effect, the Q-R algorithm? involves an 
iterative process for finding a succession of simi- 
larity transformations each of which reduces 
Amax (B) until it becomes vanishingly small. ) 


If we choose a transformation of variables of 
the form 


~1 
yeT x (91) 
to be used in equation (90) — with the intent of 
making Apax (B) small, but not necessarily zero— 
then equation (56) becomes 


Ty = ATy +g (t) (92) 


or 


~ ~1 -1 
y=T  ATy+T g(t) (93) 
Accordingly, in view of equations (87), (88) and 
(90), we may write the solution as 


y(t th) = ety (t) 


h fH as 
+f oP Di By (tea twlttD]dr (94) 
oO 


where w(t) = tg (t). In effect, then, we solve 
the original problem by transforming to the do- 
main of the new variables y before integrating. 


From equation (94) and the relation x = Ty, 2 


solution to the problem x = Ax + g(t), may be 
computed, Moreover, if T has been chosen so as 
to reduce i,,,,, (B) sufficiently, a large integration 
step h may be used without risking numerical insta- 
bility. The validity of this procedure has been 
verified in handling the following simple problem. 


Consider the RC network shown in Figure 9. 


Figure 9, RC ladder network. 


With all resistances equal and with all capacitances 
equal, the maximum integration step which could 
be taken safely using equation (88) was found to be 
about equal to the time constant, RC. In the case 
of the modified network of Figure 10, 


Figure 10. Modified RC ladder network. 


with C, split in half and the bridging resistor r 
having a very small value relative to R, the inte~ 
gration step was found to be limited by the time 
constant rCG. 


Originally, the state-variable vector used for 
both networks consisted of all the capacitor voltages. 
However, when a different combination of voltages 
was used in treating the second network — namely, 
the voltages across capacitors C;, Cz, Cg, and C5 
but the sum and difference of the voltages across 
the two half-size capacitors — it became possible 
to increase the integration step size until it was 
equal to that used in the first network, This simple 
change of variables defined a transformation which 
was sufficient to reduce significantly the maximum 
eigenvalue of B in equation (90), though not suf- 
ficient to expose all the eigenvalues explicitly. 


The lesson to be learned from the foregoing 
considerations is evidently this: one key to the 
minimum time constant problem seems to lie in 
operating on the coefficient matrix A in some 
meaningful way. But if the A matrix is left intact, 
then we must accept the limitations implicit in 
equation (84). The simplest kind of operation on 
the A matrix is to extract its diagonal terms to be 
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handled analytically as in equation (88), The next 
more difficult operation is to transform the matrix 
as in equations (90) to (94). The ultimate, of 
course, is to solve the eigenproblem completely, 


Thus, there are clear guidelines to follow in 
tackling the integration problem for nonlinear 
systems where the Jacobian matrix now plays the 
role of the A matrix. Obviously, solving the 
eigenproblem repeatedly is inappropriate, since 
the Jacobian matrix changes during the course of 
the solution, However, finding a quick and dirty" 
transformation T to be used in conjunction with 
equations (90) to (94) does seem plausible. More- 
over, it should not be necessary to update T at 
each integration step since it is likely that a given 
transformation could be used for many integrations 
before becoming degraded to the point of useless- 
ness. 


The feasibility of this approach would depend 
not only on minimizing the frequency of updating 
T but also on developing an efficient algorithm for 
computing this transformation. If the average 
computational cost per integration step can be 
kept proportional to the square of the number of 
equations, then the minimum time constant barrier 
can be breached for nonlinear problems as well 
as for linear. 


In the case of linear systems, the superiority 
of the eigenvalue approach for both a~c and tran- 
sient analysis is unmistakable, since it goes 
directly, and with reasonable computational effi- 
ciency, to the very heart of the problem. This is 
not too surprising in view of the well known facts 
that the eigenproblem is the most fundamental 
problem in matrix theory and that its solution 
provides essentially all there is of real value to 
be learned about a matrix. Fortunately, the 
advent of the digital computer, and the Q-R algo- 
rithm, now brings this powerful theoretical tool 
down to a practical working level. 


Future Developments 


It is evident from the topics discussed in this 
paper that much work needs to be done in the 


immediate future to upgrade the network analysis 
programs currently available to circuit designers, 
In addition to the incorporation of new and improved 
analytical techniques, these upgraded programs 
will need to provide facilities for statistical analy- 
sis of circuits, for storage and retrieval of device 
models, for handling transmission lines, and for 
treating much larger problems. Graphical presen- 
tation of output variables in a variety of forms will 
be essential and a much closer degree of man-~ 
machine interaction will become possible through 
the use of CRT displays and remote access to the 
computer. A good example of this kind of develap- 
ment is Katzenelson's AEDNET program. 2,40 


In another direction, there is a growing need 
to handle integrated and monolithic circuits. This 
new technology is already placing strong demands 
on our analytical, numerical, and programming 
techniques and it will inevitably force the develop- 
ment of more powerful methods for solving partial 
differential equations. This, in turn, is likely to 
stimulate the exploration and exploitation of such 
techniques as Kron's method of interconnecting 
solutions, 8-10 the use of higher dimensional net- 
work models, @! and a novel numerical method for 
treating partial differential equations recently 
described by Polozhii. * 


A step beyond the production of more powerful 
and effective network analysis programs will be 
the development of optimization techniques -- in 
which the analysis tasks will likely be relegated to 
subroutine status. Significant steps are already 
being taken in this direction. 42 Moreover, auto- 
matic network i? Cpobalied is also being carried out 
to some extent, 2 It is probable, however, that 
for a long time to come a majority of computer 
aided design will be based primarily on analysis 
and optimization techniques, with the engineer him- 
self not only playing an essential role in the "design 
loop"' but also controlling the entire process of 
design. Thus, while the computer will certainly 
take over the routine computational tasks of the 
engineer, this will enhance rather than supplant his 
creative abilities--which are really his raison d'étre. 
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